# Required packages ####
# install.packages("quantreg")
#install.packages("haven")
require(quantreg)
require(haven)

# Folders ####
directory<-Sys.getenv("US_Ineq_Repl")
MMwd<-file.path(directory, "Scripts/MM")
Datawd<-file.path(directory, "Processed")

# Import Data ####
setwd(Datawd)
Data<-read_dta('Data86-15.dta')
save(Data,file="Data86-15.Rda")
rm(Data)

# Functions ####
setwd(MMwd)
source("MMfun.R")

# Code begins ####
y0<-1986
y1<-2015
D<-GetData(size=80000,boot=F,directory=Datawd,y0=y0,y1=y1)

vars<-c("lrwage3","union","pubsect","manuf","nonwhite",
        "female","partte","married","smsa",
        "exper" , "exper2" ,
        "exp_cl" , 
        paste0("expe",1:length(unique(D$Z0$exp_cl))), 
        "years_educ" , "educ_cl" , 
        paste0("educ",1:length(unique(D$Z0$educ_cl))), 
        "edex" ,
        "region",
        paste0("reg",1:length(unique(D$Z0$region))),
        "ee_cl",
        paste0("ee",1:length(unique(D$Z0$ee_cl))),
        "ind21",
        paste0("indu",1:length(unique(D$Z0$ind21))),
        "state",
        paste0("state",1:length(unique(D$Z0$state))),
        "orgwgt","year")

## 4 regions ####
# 1 "Northeast"
# 2 "Midwest"
# 3 "South"
# 4 "West"
# ---
## 21 Industries ####
# 1 "Agriculture, Forestry, Fishing, and Hunting"
# 2 "Mining"
# 3 "Construction"
# 4 "Durable Goods Manufacturing"
# 5 "Nondurable Goods manufacturing"
# 6 "Durable Goods Wholesale"
# 7 "Retail Trade"
# 8 "Transportation and Warehousing"
# 9 "Utilities"
# 10 "Information"
# 11 "Finance and Insurance"
# 12 "Real Estate and Rental and Leasing"
# 13 "Professional, Scientific, and Technical Services"
# 14 "Management, Administrative and Support, and Waste Management Services"
# 15 "Educational Services"
# 16 "Health Care and Social Assistance"
# 17 "Arts, Entertainment, and Recreation"
# 18 "Accommodation and Food Service"
# 19 "Other Services (Except Public Administration)"
# 20 "Public Administration"
# 21 "Armed Forces"
# ---
## 50 States and DC ####
# 11 "Maine"
# 12 "New Hampshire"
# 13 "Vermont"
# 14 "Massachusetts"
# 15 "Rhode Island"
# 16 "Connecticut"
# 21 "New York"
# 22 "New Jersey"
# 23 "Pennsylvania"
# 31 "Ohio"
# 32 "Indiana"
# 33 "Illinois"
# 34 "Michigan"
# 35 "Wisconsin"
# 41 "Minnesota"
# 42 "Iowa"
# 43 "Missouri"
# 44 "North Dakota"
# 45 "South Dakota"
# 46 "Nebraska"
# 47 "Kansas"
# 51 "Delaware"
# 52 "Maryland"
# 53 "District of Columbia"
# 54 "Virginia"
# 55 "West Virginia"
# 56 "North Carolina"
# 57 "South Carolina"
# 58 "Georgia"
# 59 "Florida"
# 61 "Kentucky"
# 62 "Tennessee"
# 63 "Alabama"
# 64 "Mississippi"
# 71 "Arkansas"
# 72 "Louisiana"
# 73 "Oklahoma"
# 74 "Texas"
# 81 "Montana"
# 82 "Idaho"
# 83 "Wyoming"
# 84 "Colorado"
# 85 "New Mexico"
# 86 "Arizona"
# 87 "Utah"
# 88 "Nevada"
# 91 "Washington"
# 92 "Oregon"
# 93 "California"
# 94 "Alaska"
# 95 "Hawaii"
# ---

## Formula ####
omit<-c("lrwage3",
        "exp_cl","expe1",
        "years_educ",
        "educ_cl","educ1",
        "ee_cl","ee1",
        "region",paste0("reg",1:length(unique(D$Z0$region))),
        "ind21", "indu1",
        "state", "state1",
        "orgwgt","year")
xvars<-vars[!vars %in% omit]

f<-paste0("lrwage3 ~ ",paste(xvars,collapse =" + ")) 

formula <- as.formula(f)
h<-paste("c('(Intercept)',", gsub("\\+",",",gsub(" ","'",gsub("^.*?~","",f))),"')",sep="")
h<-eval(parse(text=h))

md<-lm(formula,data=D$Z0,weights=orgwgt)
summary(md)

omit0<-names(md$coefficients[is.na(md$coefficients)]) 
xvars0<-xvars[!xvars %in% omit0]
f0<-paste0("lrwage3 ~ ",paste(xvars0,collapse =" + ")) 
formula0 <- as.formula(f0)

qmd<-rq(formula0,data=D$Z0,tau=0.5,weights=orgwgt,method="pfn")
summary(qmd)

J<-69
taus <- 1:J/(J+1)

start_time <- Sys.time()
model0<-rq(formula0, data=D$Z0,tau=taus,weights=orgwgt)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
#time = 3.33 hours

start_time <- Sys.time()
smodel0<-summary(model0)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
#time = 6.76 hours 

md<-lm(formula,data=D$Z1,weights=orgwgt)
summary(md)

omit1<-names(md$coefficients[is.na(md$coefficients)]) 
xvars1<-xvars[!xvars %in% omit1]
f1<-paste0("lrwage3 ~ ",paste(xvars1,collapse =" + ")) 
formula1 <- as.formula(f1)

qmd<-rq(formula1,data=D$Z1,tau=0.5,weights=orgwgt,method="pfn")
summary(qmd)

start_time <- Sys.time()
model1<-rq(formula1, data=D$Z1,tau=taus,weights=orgwgt)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
#time = 1.7 hours

start_time <- Sys.time()
smodel1<-summary(model1)
end_time <- Sys.time()
di<-end_time - start_time
cat(paste0('time = ',round(di,2),' ',di%.%units),'\n')
#time = 3.43 hours

setwd(Datawd)
save.image("QR_mod.RData")
